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Abstract. Robust heteroclinic cycles are known to change stability in resonance bifurcations, 
which occur when an algebraic condition on the eigenvalues of the system is satisfied and which 
typically result in the creation or destruction of a long-period periodic orbit. Resonance bifurcations 
for heteroclinic networks are potentially more complicated because different subcycles in the network 
can undergo resonance at different parameter values, but have, until now, not been systematically 
studied. In this article we present the first investigation of resonance bifurcations in heteroclinic 
networks. Specifically, we study two heteroclinic networks in R 4 and consider the dynamics that 
occurs as various subcycles in each network change stability. The two cases arc distinguished by 
whether or not one of the equilibria in the network has real or complex contracting eigenvalues. 
We construct two-dimensional Poincare return maps and use these to investigate the dynamics of 
trajectories near the network; a complicating feature of the analysis is that at least one equilibrium 
solution in each network has a two-dimensional unstable manifold. We use the technique developed 
in |18| to keep track of all trajectories within these two-dimensional unstable manifolds. In the case 
with real eigenvalues, we show that the asymptotically stable network loses stability first when one of 
two distinguished cycles in the network goes through resonance and two or six periodic orbits appear. 
In some circumstances, asymptotically stable periodic orbits can bifurcate from the network even 
though the subcycle from which they bifurcate is never asymptotically stable. In the complex case, 
we show that an infinite number of stable and unstable periodic orbits are created at resonance, and 
these may coexist with a chaotic attractor. In both cases, we show that near to the parameter values 
where individual cycles go through resonance, the periodic orbits created in the different resonances 
do not interact, i.e., the periodic orbits created in the resonance of one cycle are not involved in 
the resonance of the other cycle. However, there is a further resonance, for which the eigenvalue 
combination is a property of the entire network, after which the periodic orbits which originated 
from the individual resonances may interact. We illustrate some of our results with a numerical 
example. 
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1. Introduction. Heteroclinic cycles and networks are flow invariant sets that 
can occur robustly in dynamical systems with symmetry, and are frequently associated 
with intermittent behaviour in such systems. Various definitions of heteroclinic cycles 
and networks have been given in the literature; for examples, see [5TfT7iri9 .20.24]. We 
use the following definitions from [18] . For a finite-dimensional system of ordinary 
differential equations (ODEs), we define: 

Definition. A heteroclinic cycle is a finite collection of equilibria {£1, . . . of 
the ODEs, together with a set of heteroclinic connections {71 (t), ■ • ■ ,ln{t)}, where 
7j(t) is a solution of the ODEs such that jj(t) — > £j as t — > — 00 and 7j(i) — > as 
t — > 00, and where £n+i = £i- 
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Definition. Let Ci,C2, . . . be a collection of two or more heteroclinic cycles. We 
say that Af = \J i Ci forms a heteroclinic network if for each pair of equilibria in the 
network, there is a sequence of heteroclinic connections joining the equilibria. That is, 
for any pair of equilibria £j , £ Af, we can find a sequence of heteroclinic connections 
{ipi (*))•••) Ipi (t)} G N an d a sequence of equilibria {£, mi , . . . , £m i+1 }sJV such that 
£ TO1 = £j, £ m , +1 = and 7 Pi is a heteroclinic connection between £ mi and £ mi+1 . 

More generally, the heteroclinic orbits in a heteroclinic cycle may connect flow 
invariant sets other than equilibria (e.g., periodic orbits or chaotic saddles) but we will 
not consider such possibilities in this article. Our definition of a heteroclinic network 
does not require that there be an infinite number of heteroclinic cycles in a network, 
but in the networks we consider, (at least) one of the equilibria in the network has a 
two-dimensional unstable manifold and associated with this is an infinite number of 
heteroclinic connections between that equilibrium and another. We only consider the 
case that the set of equilibria in the network is finite. 

Methods for determining the stability properties of an isolated heteroclinic cycle 
involving equilibria or periodic orbits are well-established [TTl[T9j[21 23,26,22], and 
their implementation is generally straightforward, at least in principle, because there 
is only one route around the cycle. In the most widely studied examples, all equi- 
libria have one-dimensional unstable manifolds, and within these manifolds, the next 
equilibrium point in the cycle is a sink. One way a heteroclinic cycle can lose stability 
is in a resonance bifurcation. A resonance bifurcation is a global phenomenon, which 
occurs when an algebraic condition on the eigenvalues of the equilibria in the cycle is 
satisfied. Gcnerically, resonance bifurcations are accompanied by the birth or death of 
a long-period periodic orbit. If the bifurcation occurs supercritically, then in the sim- 
plest case, the bifurcating periodic orbit is asymptotically stable and the heteroclinic 
cycle changes from being asymptotically stable to having a basin of attraction with 
measure zero. Resonance bifurcations from asymptotically stable heteroclinic cycles 
have been extensively studied; see [T21I22, 26, 27]. for cases in which all eigenvalues 
are real, and [25] for a case with complex eigenvalues. Much less is known about 
resonance bifurcations of non-asymptotically stable cycles. 

Stability of robust heteroclinic networks is less well understood. Some results are 
known (e.g., [3ll5 HTfJl[LH[TMT51l20H24) ) but these are, in general, partial results and 
confined to specific examples. One source of difficulty is that there may be many 
different routes by which an orbit can traverse a heteroclinic network, and keeping 
track of all possibilities in the stability calculations can be challenging, particularly 
when one or more of the equilibria in the network has a two-dimensional unstable 
manifold. When this occurs, trajectories may go straight to an equilibrium point 
that is a sink within the unstable manifold, or may visit a saddle equilibrium point 
before moving on to the sink. A full analysis needs to account for all possibilities. 
In [T5] , we showed that it is possible to do this and so to establish relatively complete 
stability results for a specific class of problems in which all cycles in the network share 
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a common hetcroclinic connection, despite there being several equilibria with two- 
dimensional unstable manifolds. In this case, we were able to derive conditions that 
determine the attractivity properties of the network. These conditions are network 
analogues of stability conditions for single hetcroclinic cycles, and involve inequalities 
on combinations of the eigenvalues of the equilibria. By analogy with resonance 
bifurcations of heteroclinic cycles, we call the transition that occurs when one or 
more of the inequalities is reversed a resonance of a heteroclinic network. 

In |18j it was noted that complicated dynamics could be associated with resonance 
in the network studied. Our aim in this article is to complete this analysis, and extend 
it to a closely related heteroclinic network (which is the same as that studied in (TT] . 
although in that article, no attempt was made to keep track of all trajectories in the 
two-dimensional unstable manifolds). We will then investigate resonance bifurcations 
in both networks in detail. We believe this is the first article to analyze network 
resonance in a systematic way. 

Both networks have the basic structure shown schematically in figure fTTTl Specif- 
ically, each network consists of six equilibria, which we call A, B, X, Y, P and Q, and 
their symmetric copies, —A, —B, —X, —Y, —P, and — Q, along with a collection of 
hetcroclinic connections between equilibria. The equilibria A and B are connected by 
a single (one-dimensional) heteroclinic connection from A to B. Equilibrium B has a 
two-dimensional unstable manifold associated with two different positive real eigen- 
values, and there is a continuum of heteroclinic orbits lying within this manifold and 
connecting B to X, Y, P, Q, and their symmetric copies: there is a single connection 
from B to P and from B to Q, but an uncountable family of connections from B to X 
and from B to Y. The equilibria X and Y have one-dimensional unstable manifolds 
which are heteroclinic connections to A and —A. P and Q have two-dimensional un- 
stable manifolds consisting of single heteroclinic connections to X and Y (and their 
symmetric copies) and continua of heteroclinic connections to A and —A. 

The feature that distinguishes our two networks from one another is whether or 
not the Jacobian matrix of the flow evaluated at A has complex eigenvalues. In Case I, 
there are only real eigenvalues at A, while in Case II, A has a complex conjugate pair 
of eigenvalues with negative real part. Further details about the networks are given 
in Section [2] 

We analyse the networks by deriving local and global maps that approximate the 
dynamics near and between the different equilibria in the network. This analysis is 
complicated by the fact that, for reasons explained in detail later, it is not always 
possible to write these maps explicitly. However, under certain approximations and 
assumptions about the dynamics near the networks, we are able to compose the maps; 
these approximations and assumptions mildly restrict the validity of our results. This 
then gives us information about the dynamics of all possible trajectories as they 
traverse the network and return close to where they started. The derivations of the 
maps, approximations and compositions are contained in sections [3] and 01 
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Using the return maps, we are then, in section [5J able to determine existence 
criteria for fixed points of the maps, which correspond to periodic orbits in the original 
flow. These periodic orbits appear when resonance conditions for the network are 
broken. In the case of an asymptotically stable network losing stability, we find that 
the first conditions to be violated are those associated with one or the other of the 
subcycles within the network, that is, the conditions on the eigenvalues arc the same 
as for a single cycle. In the network with real eigenvalues, either two or six periodic 
orbits appear at this initial resonance (including all symmetric copies). We also show 
that an asymptotically stable periodic orbit can bifurcate from a non-asymptotically 
stable heteroclinic cycle in this network. In the network with complex eigenvalues, 
we find that infinitely many periodic orbits appear at resonance. For both networks, 
if we remain in parameter space close to the point where the resonances of individual 
subcycles occur (we consider the eigenvalues of the equilibria as parameters), then the 
periodic orbits arising from the bifurcations of the subcycles do not interact, i.e., the 
periodic orbits created in the resonance of one cycle are not involved in the resonance 
of the other cycle. However, we find that there is a further resonance, for which the 
eigenvalue combination is a property of the entire network, after which the periodic 
orbits which originated from the individual resonances may interact, for instance when 
orbits arising from different resonances come together in saddle-node bifurcations. 

In addition to bifurcating periodic orbits, we also find that a chaotic attractor 
may be created at a resonance bifurcation of the network with complex eigenvalues. 
This is detailed in section 15.2.41 Section 15.2.51 contains a numerical example showing 
both periodic orbits and a chaotic attractor. 

In section |6] we look at resonance bifurcations of an isolated heteroclinic cycle 
with complex eigenvalues. When this cycle goes through resonance, infinitely many 
periodic orbits appear, in a similar manner to that seen within the network with 
complex eigenvalues. The analysis of this cycle allows us to conjecture which features 
of the dynamics of our Case II network arise from the existence of complex eigenvalues 
and which are a result of the network structure. 

Section [7] concludes with discussion and avenues for further work. 

2. The heteroclinic networks. We consider a system of ordinary differential 
equations, x = f(x), where x = (x±,X2, x 3 , y 3 ) € R 4 and f : IR 4 —> R 4 is a C 1 
vector-valued function. For both networks we consider, we assume this system has 
the following cquivariancc properties: 



«i(f(x)) = f(Ki(x)), 



i = 1,2, 



where 



«i : (x 1 ,X2,x s ,y 3 ) -> (-xi,x 2 ,x 3 ,y 3 ), 
k 2 : (x 1 ,x 2 ,x 3 ,y 3 ) -> {x 1 ,-x 2 ,x 3 ,y 3 ). 



(2.1) 
(2.2) 
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Fig. 1.1. Schematic diagram showing the basic structure of the two heteroclinic networks stud- 
ied. For clarity, the equilibrium Q is not shown, but this equilibrium has a similar role to equilibrium 
P except that there is a heteroclinic connection from Q to —X instead of Q to X, where —X is a 
symmetric copy of X. The remaining (conjugate) parts of the network are obtained under the ac- 
tion of the symmetry groups described in Section^ Thin curves represent single (one- dimensional) 
heteroclinic connections and bold curves represent a two-dimensional family of connections between 
the relevant equilibria. The double arrowhead on the connection from B to X indicates that expan- 
sion near B in the direction of this connection is stronger than expansion in the direction of the 
connection from B to Y . 

In Case I, we further assume that the system is equivariant with respect to the sym- 
metries 



while in Case II we assume that the system is also equivariant with respect to the 
symmetry 



Note that the symmetries K\, k 2 , k-x and n y are those used in the network in j!7j 
while the symmetries Ki, k 2 and n 3 are those used in the network in [18]; imposing 
the assumptions listed below ensures that Case I is precisely the network from |17) 
and Case II is the network from [T5] . 

The equivariance properties of the networks cause the existence of dynamically 
invariant subspaces in which robust saddle-sink heteroclinic connections can occur. 
We make the following further assumptions about the dynamics in these subspaces, 
as illustrated in figure |2~T1 

• Al: There exist symmetry-related pairs of equilibria ±A and ±B on the X\ 
and X2 coordinate axes, respectively. Within the invariant plane x 3 = y 3 = 0, 
A is a saddle and B is a sink and there is a heteroclinic connection from A 
to B. See HaureWlVa). 

• A2: There exist symmetry-related pairs of equilibria ±X , ±Y , ±P and ±Q 
in the invariant plane x\ = x 2 =0. Within this subspace, ±X and ±Y are 
sinks, while ±P and ±Q are saddles. The eight equilibria together with the 



n x : (x 1: x 2 ,x 3 ,y 3 ) -> {xi,x 2 ,-x 3 ,y 3 ), 
k v : (x 1 ,x 2 ,x 3 ,y 3 ) (x 1 ,x 2 ,x 3 ,-y 3 ). 



(2.3) 
(2.4) 



«3 : (xi,x 27 x 3 ,y 3 ) (x!,x 2 ,~x 3 ,-y 3 ). 



(2.5) 



5 




Fig. 2.1. Dynamics within invariant subspaces of the two networks being considered. For 
clarity, only part of the relevant subspaces are shown in panels (a), (c), (d) and (e), with the 
dynamics in the omitted parts being obtained by applying the symmetries, (a) The invariant plane 
X Z = Vs = 0, showing the heteroclinic connection from A to B. (b) The invariant plane x\ = £2 = 0, 
showing the invariant circle C and the equilibria iX, ±Y, ±P and ±Q that lie on C. (c) The 
subspace x\ = showing part of the two-dimensional unstable manifold of B and part of the circle C 
in the (£3,1/3) plane. The equilibria ±X and ±Y are shown to lie on the coordinate axes, with the 
eigenvectors of the corresponding linearised flow at B aligned with the axes. In the Case I network, 
this situation is forced by the symmetries k x and n y . In Case II, ±X and ±Y are chosen, for 
convenience, to lie on the axes, but no assumption is made about the alignment of the eigenvectors, 
(d) The subspace X2 = for Case I, showing connections from X, Y and P to A. The unstable 
manifold of Q (not shown) behaves similarly to the unstable manifold of P. The connection from X 
to A (resp. Y to A) lies in the invariant plane X2 = y?, = (resp. X2 = £3 = 0). (e) The subspace 
X2 = for Case II, showing spiralling of the unstable manifolds of X, Y and P into A. The unstable 
manifold of Q (not shown) behaves similarly. In each subspace, the flow is strongly contracting in 
the radial direction. 



heteroclinic connections between them make up an invariant curve C , which 
is topologically a circle. We hereafter refer to C as a circle, and we assume 
that C can be parametrised by the angle 63, the polar angle in the (23,2/3)- 
plane. Note that the intersections of the stable manifolds of ±P and ±Q 
with the invariant plane form the boundaries between the basins of attraction 
of ±X and ±Y in the invariant plane. Only a small part of each intersection 
is shown in figure l2TT b). to avoid giving a misleading impression about the 
dynamics near the origin of the (2:3, i/3)-plane, but each intersection curve in 
fact extends to the origin of the subspace. In Case I, the X3 and 1/3 axes are 
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invariant and coincide with orbits of the system, but this is not necessarily 
so in Case II. 

• A3: Within the invariant subspace x\ = 0, there exist two-dimensional mani- 
folds of saddle-sink connections from B to ±X and ±Y ffiaure [2AV c)). There 
are also one- dimensional (saddle-saddle or saddle-sink) heteroclinic connec- 
tions from B to ±P and ±Q and from ±P and ±Q to ±X and ±Y , as 
shown in Haure \2. lY c). The unstable manifold of B is two-dimensional, and 
the stable manifolds of ±X and ±Y are each three-dimensional within the 
subspace. In Case I, there is a connection from B to X (resp. from B to Y) 
in the subspace X\ =2/3 = (resp. the subspace X\ = X3 = 0). 

• A4: Within the invariant subspace x-i = 0, there exists a two-dimensional 
manifold of saddle-sink connections from C to A. Within this manifold, A is 
either a stable node (Case I) or a stable focus (Case II). A similar manifold 
connects the equilibria on C to —A. Apart from the heteroclinic connections 
from ±P and ±Q to ±X and ±Y, the unstable manifolds of ±P and ±Q are 
contained in the stable manifolds of A and — A. There are no equilibria other 
than the origin and those mentioned above lying in the subspace x-i = 0. See 
fiaure [KlV d) and (e). 

• A5: Equilibrium B has real eigenvalues corresponding to dynamics in its 
unstable manifold, and these eigenvalues are unequal. We do not consider the 
case where B has complex eigenvalues. 

Assumptions Al A5 ensure the existence of the two heteroclinic networks con- 
sidered in this article. The symmetries K\ and K2 ensure that X\ and X2 cannot change 
sign along a trajectory, so we consider x\ > and xi > only. Similarly, in Case I, 
the symmetries k x and n y ensure that X3 and 2/3 cannot change sign along a trajectory, 
so in this case we can consider X3 > and 2/3 > only. 

To simplify our analysis, we make several further assumptions. The first part of 
A7 and assumptions A8 and A9 are automatically satisfied for Case I, but we extend 
them to Case II as well. A6 is a genericity assumption and A8 is not restrictive. 
Either A7 or A9 can always be satisfied; we restrict the dynamics by assuming both 
are true. A10 is a restrictive assumption. 

• A6: At each equilibrium, no two of the eigenvalues of the linearisation are 
equal. 

• A7: The two expanding eigenvectors at B lie in the X3 and 2/3 directions. 
Without loss of generality we assume that the eigenvalue with eigenvector 
pointing in the X3 direction is larger than that corresponding to the 2/3 direc- 
tion. 

• A8: The linearisation around A is in Jordan form. 

• A9: The equilibria ±X and ±Y are, respectively, on the X3 and 2/3 coordinate 
axes. 

• A10: At A (resp. B) the strong stable direction lies along the coordinate axis 
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Xi (resp. x-i). At each of ±X , ±Y , ±P and ±Q, the strong stable direction 
lies in the (x^,y^) plane and is transverse to C (which is an invariant circle, 
by A2). 

We can therefore summarise the different networks we study as follows. 

• In both cases, the overall network is A — > B — > C — > A, where, within C, 
trajectories can visit any of ±A, ±F, ±P and ±Q, although only in certain 
orders as indicated in figures 11.11 and 12.11 All cycles in the network contain 
either three or four equilibria. 

• Case I is equivariant under the symmetries m, K2, k x and n y , and the lin- 
earisation of the vector field at each equilibrium has only real eigenvalues. 

• Case II is equivariant under the symmetries k 1; k 2 an d k 3 , and the lineari- 
sation at equilibrium A has a pair of complex conjugate eigenvalues with 
negative real part. 

3. Maps for the dynamics near the heteroclinic networks. We follow 
the standard procedure for modelling the dynamics near a heteroclinic network, i.e., 
we construct return maps defined on various cross-sections in R 4 and analyze the 
dynamics of these maps. Cross-sections transverse to the connection from A to B 
are of special interest, since all trajectories lying near one of our networks must pass 
through such a cross-section and so maps defined on such a cross-section contain 
information about the asymptotic stability of the network as a whole. However, in 
our investigation of resonance bifurcations, it will be important to consider situations 
in which the network has more subtle stability properties, in which case we will be 
interested in return maps defined on cross-sections transverse to other heteroclinic 
connections. 

In Section 13.11 we give details of the coordinates, cross-sections, and local maps 
(valid near equilibria) we use in construction of the return maps. Apart from the local 
map near A in Case I, these are the same as the maps found in [18]. In Section l3~2l 
we derive global maps (valid near heteroclinic connections between equilibria); these 
are the same as the global maps found in |18| apart from some additional constraints 
needed for Case I. The local and global maps we define are consistent with those 
used in [17], but have a more general form (and use different notation), since here 
the maps are designed to capture the behaviour near the whole heteroclinic network, 
whereas in [T7] the analysis focussed on two distinguished cycles (called the ^3-cycle 
and the ^4-cycle in [17] , corresponding to the heteroclinic cycles through X and Y in 
the notation of this article). 

In principle, the local and global maps can be composed in an appropriate order 
to obtain return maps modelling the dynamics near our networks. However, because 
we wish our maps to keep track of a continuum of heteroclinic cycles in network, it 
turns out that we are unable to derive explicit forms for some of the local maps and 
hence for the return map as a whole. However, we are able to obtain approximations 
of the maps for particular ranges of the coordinates in our maps, and this is sufficient 
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for us to be able to extract results about resonance. 



3.1. Coordinates, cross-sections, and local maps. Near A and B, we define 
local coordinates that place the equilibrium at the origin. We write Xj or j/j if the 
local coordinate is the same as the corresponding global coordinate, and use Ui for 
the local coordinate otherwise. We use polar coordinates when it is more convenient: 
{x 3 .,y 3 ) becomes (73,03), where x 3 = r 3 cos9 3 and 2/3 = r 3 sin0 3 , and U3 measures 
the distance within the 2:3-2/3 subspace from the invariant circle C. Assumptions A7 
and A8 guarantee that the coordinate axes are aligned with the eigenvectors of the 
relevant linearised system. 

Near A, the linearised flow in Case I is given by: 

ui = -r A ui, x 2 = e A x 2 , ±3 = -c Ax x 3 , y 3 = -c Ay y 3 , (3.1) 

where Ta, sa, cax and ca v are positive constants. The letters e, c and r in these 
constants refer to the expanding, contracting and radial directions, as defined by [21] . 
In Case II, the linearised flow near A is given by: 

Ui = -r A Ui, x 2 = e A x 2 , ± 3 = -c A x 3 - uy 3 , y 3 = ujx 3 - c A y 3 , (3.2) 

where ta, e^, ca and 10 are positive constants. In polar coordinates, the x 3 and 2/3 
equations give r 3 = —cat?, and #3 = cu. 
Cross-sections near A are defined as: 

= {(u 1 ,x 2 ,r 3 ,e 3 )\\u 1 \<h,0<x 2 <h,r 3 = h,0<e 3 <27r}, 
H^ ut ee {(ui,X2,r 3 ,6 3 )\\u 1 \<h,x 2 = h,0<r 3 <h,0<6 3 <27r}. 

Here < h <C 1 is a parameter small enough that the cross-sections lie within the 
region of approximate linear flow near A (and similarly near B and C, as required 
below) . 

In Case I, the flow near A induces a map cj>A.r ■ HJ4 —> H^ ut , which is obtained 
to lowest order by integrating equations (13. 1|) : 



4>A,r{ui,X2,h,9 3 ) 



( —J , h, h \^cos 9 3 + sin 9 3 ^ 

X 2 \ S A V -$A 



tan" 1 tan 9 3 (-±) , (3.4) 



where S Ax = ^f^ and 8 A y = -^ JL - In Case II, the corresponding local map, obtained 



to lowest order by integrating equations (|3.2j) . is 



(«i ,X2,h,0 a )= U (^) eA ,h,h (^) A ,9 3 - 2- log (^) J , (3.5) 



where S A = — ■ 



Near B, the linearised flow is: 



±i = -c B x 1 , u 2 = -r B u 2 , ±3 = e Bx x 3l y 3 = e By y 3 , (3.6) 

where r B , e Bx , e By , c B are positive constants. From A7, we have e Bx > e By . Cross- 
sections near B are defined as: 

H£ = {(x u u 2 ,r 3 ,9 3 )\x 1 =h,\u 2 \<h,0<r 3 <h,0<9 3 <2ir}, 
H^ ut = {{x l ,U2,r 3 ,9 3 )\Q<x 1 <h,\u 2 \<h,r 3 = h,Q<6 3 <2-K}, 

and the flow induces a map <p B : 1V B — > H™', which is obtained to lowest order 
by integrating equations (|3.6I) . The map cannot be written down explicitly, but is 
computed as follows. First, the x 3 and y 3 equations are solved: 

x 3 {t) = r 3 {0)cos8 3 {Q)e eB *\ y 3 (t) = r 3 (0) sin# 3 (0) e eB »*, 

where r 3 (0) and 3 (0) are the initial values of the radial coordinates (i.e., on H 1 ^). 
The trajectory crosses H| ut when r 3 (t) = h, so the transit time T B is found by solving 
the equation 

h \ 2 

— - ) = cos 2 9 3 (0) e 2eB * TB + sin 2 3 (O) e 2eB * TB (3.8) 

73(0)7 

for T B in terms of r 3 (0) and 9 3 (0). This yields the local map (f> B : Hg H™*: 

0B(/i,U2,r 3 ,6l3) = (he~ CBTB ,u 2 e- rBTB ,/i,tan _1 (tan6» 3 e (e ^- e ^ )TB )) . (3.9) 

For later convenience, we define <5ro- = - £a - and S Bv = 

Near the circle C we would like a local map that captures the dynamics of all 
orbits that pass near C. Linearization of the flow near the equilibria on C alone will 
be insufficient for our purposes. Instead, we use the technique described in j!8j and 
summarised below to construct a map. Specifically, we assumed in A2 that C can 
be parameterised by the angle 8 3 . The rate of relaxation onto C is controlled by the 
^-dependent quantity — rc(0 3 ). The assumption of strong contraction in the radial 
(r 3 ) direction (A10) means that the dynamics on C of 9 3 can be captured by an 
equation of the form 9 3 = g(9 3 ), where g is a nonlinear function with g(0) = g(^) = 
ff( 7r ) = ff(nr) = (this last statement follows from assumption A9 which stipulates 
that ±X and ±Y lie on the coordinate axes). There will be further zeroes of g at the 
values of 9 3 corresponding to ±P and ±Q. These considerations mean we can model 
the flow near C by: 

ii = e c {9 3 )xi, x 2 = -c c (9 3 )x 2l u 3 = -r c {9 3 )u 3l 9 3 = g(9 3 ), (3.10) 

where re, ec and cc are positive functions of 9 3 . 
Cross-sections near C are defined as: 

HjJJ = {(xi > x a ,tt3.03)|O<xi<Ma = Mu3|<M<03<27r} ) 
H£ ut = {(x 1 ,x a ,tt3.03)|a:i = M<S2<Mu3|<fc,O<03<27r}. 
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(3.11) 



The local flow near C induces a map 4>c '■ H£ — > H™*. We cannot write down 
the map explicitly, but it is computed as follows. First, the 8 3 equation is solved using 
an initial condition #3(0), yielding 3 (t). Then the x\ and X2 equations are solved: 

Xl {t) = x 1 (0)exp Qf ec(6 3 (t'))dt'^ , a: 2 (t) = hcxp (- J c c (e 3 (t')) dt'^j . 

The trajectory crosses H™' when x\(i) = h, so the transit time Tc can be found by 
solving 

\c(8 3 (t'))dt' = - log (3.12) 

for Tc in terms of the initial values a;i(0) and #3(0) on Then the local map 

<t>c : Hg -> H° ut is given by 



4> c {x u h,u 3 ,6 3 ) = U^expf-jT C c c (0 3 (i'))^'j , W 3(Tc),03(T c )j , (3.13) 
where U3(Xb) = 1*3 exp 

(- J Tc rc(6 3 (t')) df) . For later convenience, we define 5cx 
and 5cy to be the ratio ^(83) evaluated at the points X and Y respectively. 

As noted above, neither of the maps <pB and <pc can be written down explicitly. 
In the case of <pB, this is because we cannot write down an explicit solution of (|3.8I) 
for the transit time Tb- In the case of <f>c, the nonlinear evolution of #3 on C is not 
known explicitly. In section 14.11 we make assumptions about the flow near C and are 
then able to make approximations to the local maps in order to compute stability and 
bifurcation properties of the network. 

3.2. Global maps. To construct global maps that approximate the dynam- 
ics near heteroclinic connections of the networks, we linearise the dynamics about the 
unstable manifold leaving each of A, B and C. In doing so, we allow for the fact that 
the unstable manifold of A is one-dimensional, but the unstable manifolds of B and C 
are two-dimensional. The different equivariance properties of the vector fields for our 
different networks result in different constraints on the global maps for Case I and II. 

The heteroclinic connection from A to B intersects H™* at {u\,X2,x 3 , y 3 ) = 
(0, h, 0, 0), and intersects H'JJ at (xi, U2, x 3 , y 3 ) = (h, es, 0, 0), for a small constant es- 
Without loss of generality, we assume that es 7^ 0. Here and below, the e parameters 
give the value of the local radial coordinate at the intersection of the heteroclinic 
connection with the incoming section. These turn out to play no role at leading order, 
which is consistent with results about radial eigenvalues for heteroclinic cycles [2Tj . 

Generically, the dynamics near the heteroclinic connection will be (to lowest order, 
and in cartesian coordinates) an affine linear transformation. In polar coordinates, 
this yields, at leading order: 



*AB(ui,h,r a ,9a) = (h,e B ,D B (0 3 )r 3 ,6 B (e 3 )), 
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(3.14) 



where Dg(9 3 ) is an order-one function of 9 3 and 9g(9 3 ) is an order-one function of #3. 
Invariance of the map under the symmetry n 3 (for Case II) has the same effect on the 
form of the map as invariance under k x and K y (Case I), i.e., it ensures that there is 
no constant term or linear dependence on Ui in the ^-component. Thus, the form of 
^ ab given above is valid for both the heteroclinic networks we consider. However, 
in Case I, the invariance of the x 3 and 2/3 coordinate planes requires some additional 
constraints on the function 9b{9 3 ). Specifically, in Case I, #b(0) = 0, #b(-§) = f , 
^b(tt) = 7T and 9b(^~) — In both cases, the overall effect of the map <P ab is to 
multiply the small variable r 3 by an order-one function of #3, and to map the outgoing 
angle #3 to an incoming angle 9b{9 3 ). 

The two-dimensional unstable manifold of B intersects H™' at (x±, U2, r 3 , #3) = 
(0,0, M3) forO < 3 < 2tt, and intersects Hg at (x 1 ,x 2 ,u 3 ,9 3 ) = (0, h, e c (6» 3 ), 9 C (9 3 ))§ 
where ec is a small function of 9 3 and 9c is an order-one function of #3. To leading 
order in x\ and 112, we find: 

*bc(xi, u 2 , h, 9 3 ) = {D c (9 3 )xi, h, e c (6 3 ), 9 C {9 3 )) , (3.15) 

where Dc(9 3 ) is an order-one function of 9 3 . As for the map *&ab, in Case I there 
are additional constraints on the function 9c due to the invariance of the coordinate 
axes. Specifically, in Case I, 9 C (0) = 0, c (f) = f , M 71 ") = tt and 9 C {^-) = \- 
In both cases, we assume without loss of generality that tc{9 3 ) 7^ for any 9 3 . The 
function ec(9 3 ) plays a similar role to the constant es in (|3.14p . except that it takes 
on a different value for each heteroclinic connection and so is a function of 9 3 . In both 
cases, the overall effect of (|3.15[) is to multiply the small variable x\ by an order-one 
function of 6*3, and to map the outgoing angle #3 to an incoming angle 9c{9 3 ). 

The unstable manifold of C is two-dimensional; it intersects H™' along the 
curve (xi,X2, u 3 ,9 3 ) = (h,0,0,9 3 ), where < 9 3 < 2ir, and it intersects H™ at 
(u\, X2, r 3 , 9 3 ) = (eA(9 3 ), 0, h, 9a(9 3 )), where ca is a small function of 6*3 and 9a is an 
order-one function of 9 3 . For small X2 and u 3 , we have: 

* C A(h, x 2 ,u 3 , 8 3 ) = (e A (6 3 ),D A {9 3 )x 2 ,h, 9 A {6 3 )) , (3.16) 

where Da{9 3 ) is an order-one function of #3. In Case I, invariance of the coordinate 
axes means that 9a(0) = 0, 0a(§) = § , ^a(tt) = tt and 9a{^-) — In both cases, 
the overall effect of ^ ca is to multiply the small variable X2 by an order-one function 
of #3, and to map the outgoing angle 9 3 to an incoming angle 9a- 

4. Preliminary analysis of maps. In order to make further progress, it is 
necessary to introduce some approximations and simplifications to the local maps. 

In section |4~T1 we construct approximations to the local maps near A and B valid 
close to the X and Y directions. We also assume a simple form for the dynamics 
near C; we believe that this simplification will not qualitatively change our results. 
Throughout this section, we set h — 1 without loss of generality; this is equivalent to 
rcscaling the local coordinates introduced in the previous section. 
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Once the approximations arc made, we are then (in sections 14.21 and 14 .4)) able 
to compose the maps and compute a quantity we call 5(93) which gives the rate of 
contraction or expansion of trajectories near the network, as a function of the coor- 
dinate 63. This quantity plays a similar role to the ratio of contracting to expanding 
eigenvalues used to determine stability of some heteroclinic cycles. However, because 
we are working with a network, the ratio is dependent on the particular route taken 
around the network. As part of these calculations, we find it useful to define: 

5x = SaxSbxScx, 5 y = SAySsyScY 

where in Case II, 5ax = 5 Ay = 5 a- 

4.1. Approximate local maps. First we look at the map for the dynamics 
near A in Case I, 4>a '■ (v>i, %2> 1> $3) ~^ ( u ij 1 ; r 3i 9 3 ). ln the following, the notation 
9 Ain (resp. 6* Aout ) refers to the value of 9 3 on H™ (resp. H™*), while X2 (resp. r 3 ) 
represents the value of the second (resp. third) coordinate on H™ (resp. H™*). Then, 
from (13.41). we have 



r 3 = [cos 2 9 A ™ xf A * + sin 2 9 A ™ x? A » ) (4.1) 
and 

tan0 A °>" = tan 6^'" x S 2 Ay ~ Sji % (4.2) 

where r 3 and X2 are both small. When 9 Ain = (resp. ^), we have logr 3 = 8ax logX2 
(resp. logr 3 = 5 Ay logx 2 )- 

Expression (|4.1[) can be rewritten as 

r 3 = |cos# A -| xfr' (l+tan 2 e A ^xf SAv - SAa) y /2 

so we have 

logr 3 log Icos^l 1 log (l + tan2 Ai * xf 6 *"-'^ 



1/2 



5 Ax + 



logx 2 log £2 2 logx 2 

Note that the term inside the logarithm may be large or small. We further approxi- 
mate this later as appropriate. 

In the case of complex eigenvalues at A, the local map (|3.5[) gives: 

bgr3 5 A , e A ^=e A *-^lo g X2. (4.3) 



log X2 eA 

Approximating the local map at B is complicated by the need to solve (|3.8[) for 
the transit time Tg. At B, we have by assumption A7 that ebx > ^By, and so 
5bx < 5b v - Let 9 Bia (resp. 9 Bmt ) be the value of 9 3 on Hg (resp. H™') and denote 

13 



by T3 (resp. Xi) the value of the third (resp. first) coordinate on (resp. H™ 4 ). We 
can then rewrite (13.81) as: 



cos 



2 9 B ln £ 2e Bx T B ^ + tan 2 Q B in e 2(e By -e Bx )T E 



where r 3 is small. 

As long as 9 Bia is not too close to -| or the second term in the brackets 
is small compared to the first; we drop this term and solve for Tb, finding Tb = 
— logr3| cos6> Bin |. The term that was dropped is small (with this value of Tb) so 
long as | cot 9 Bin | > f B , where 



< 1 



(5 B y > S B x)- 



(4.4) 



When | cot8 BiI1 \ <C eB (i.e., 9 Bia is close to \ or 4^), we cannot drop the second term 
but instead approximately solve (|3.8|) . finding 



T h 



— log ( r 3 | sin C Bi 

e By 



1 

1 + - cot 



2 a B 



,2 1- 



From these expressions, we can use 
after cf>B- 



r% sin f 

to find the exit values of xi and #3 



logxi 
logr 3 



log I cos 6> S H 



>Bx 



>Bx- 



log r 3 



I cot 6» S H > 9 tB , 



5b v + 5 



cot 



2 aB in 



2(1- 



(4.5) 



By 



and 



tan ( 



tanr 



(T 3 I cos 

5 By _, 



I l + i cot 2 J 



log 7-3 



2(1- 



COt^ 



-1 



|cot0 Bi »| 

I cot e B -A 



> 



(4.6) 

where X\ and r 3 are both small. 

There are three obstacles to estimating the local map near C: the #3 dynamics is 
given by #3 = 5(6*3), where (7(6*3) is unknown, and ec(6*3) and cc{0^) are unknown. In 
order to make progress, we take simple forms for 5(6*3), ec(6*3) and 00(83) that allow 
us to solve for 6*3 (i) and to compute the required integrals. We believe that these 
simplifications will not qualitatively change our results. 

In the following we let 6 Cin (resp. 6* Cout ) be the value of 6* 3 on Hg (resp. H™ 1 ) 
and denote by x\ (resp. X2) the value of the first (resp. second) coordinate on H™ 
(resp. U° c ut ). 

We first assume that ec does not depend on 6*3. This allows us to calculate the 
transit time from H™ to H™*: 

1 



Tr 



ce- 
ll 



logxi 



We then assume that g takes a very simple form, i.e., we choose g{9 3 ) = — j sin(4# 3 ), 
with A > 0. Then X and Y are at 8 3 = and |, and P is at 9 3 = |. With this form 
for g we can solve 9 3 = g(9 3 ), and find 

tan 20 3 (i) = tan26» Cln e" At , 

taking 6> 3 (0) = (9 C "\ With 3 (T C ) = 6» C °"S we find 



tan26> 



Cout . 



tan 29 



(4.7) 



It would be tempting to assume also that cc does not depend on #3; however, this 
turns out to be too restrictive. Instead, we write 



, a v c C x + ccy . ccx - CCY nQ 
cc{9?,) = H cos26» 3 ; 



(4.8) 



this ensures cc(0) = cc(tt) = ccx and cc(f ) = cc(^) = ccy ■ With this, the exit 
value of X2 is exp (^— Jq° cc{9 3 (t')) d£ \ . From above, we know tan26>3(t) explicitly, 
so cos(26>3) = ±(1 + tan 2 (20 Cin )e _2A *) -1 / 2 , where we take the positive square root if 
< 9 CirL < f and the negative square root if \ < 9 Ci " < f . Note that 



1 



Vl + K 2 e- 



2Xt 



-At 



2A 



log 



Vl + K 2 e- 2Xt - 1 



^l + K 2 e- 2Xt + 1 



and so we find 



log £2 Sex + Scy . c C x - c C Y , 

-±^T~, log 



' ( V 1 + tan 2 2(9 C >» x{ c - l^j (| sec 26> c '" | + l) \ 



logzi 



4 A log x\ 



1 + tan 2 26» c - x{° + 1 (| sec 20 c ™ \ - 1 



7 



where X2 and x\ are both small. As before, the plus sign is taken if < # Cin < f and 
the minus sign is taken if ^ < 9 l " < ^. 

2A 

If we are away from 9 CiTL = j, such that tan 2 29 Cia x{° <C 1, then we can approx- 
imate the function above as: 



logx 2 
logxi 



' fcx-^^log(l-tan 2 ^."), O<0°*<* 
\{5cx + S C y), 



5C ta _ 7T 
7 4 ' 



(4.9) 



jCY 



+ C -%&7 10 § C 1 ~ cot2 ° Gln ) ' f < 0Ci " ^ I' 



where x 2 and Xi are both small, and the bounds near ^ are taken to mean that 

A 

t&n29 c ™x{ c < 1. 
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4.2. Composing the maps: Case I. In this subsection we consider the return 
maps for Case I. In section 14.2.11 we compose the maps starting on each of 
H'g and H™, and for each return map, we focus on the 6*3 component. We argue 
that in the parameter regimes of interest, the return maps give the same dynamics 
regardless of which section we start on. Thus in section 14.2.21 where we consider the 
other component of the return map, we need only consider the return map starting 
on H™. Note that away from resonance when the network as a whole is attracting, 
this is not the case — in order to fully describe the dynamics of trajectories near 
the network, the composition of the maps must be considered starting on all three 
Poincare sections. This observation was made in [17] and more details can be found 
in that article. A second example of this behaviour was also seen in [24] for a more 
complicated heteroclinic network. 

4.2.1. 6*3 component. As in the previous section, we denote by 9 Aia (resp. 6 Aaut ) 
the value of #3 on HJ^ (resp. H^ ut ), and by 9 Airl the value of 6*3 after one application 
of the return map from to itself; 9 Ain will typically depend on 9 Ain and x 2 . The 
symbols 9 B ' ln and 9 Cin are defined in an analogous way on the cross-sections H 1 ^ and 
H™. Without introducing ambiguity, we also write X2 for the value of X2 on H™, 
for the value of r% on and X\ for the value of x\ on H™. 

We wish to compute the derivative of 9 Ain with respect to 9 AiD at two special values 
of 6*3 , those corresponding to the invariant subspaces containing the heteroclinic cycles 
through X and Y, and similarly for derivatives of 9 in and 9 Cin . We can compute these 
derivatives without computing the entire return map, and doing so greatly simplifies 
the computation (which we give below) of the return map for general values of #3. 
Simple calculations following from section [44] give 

I x ocy-ock A in = 



d0 A - 




7 2 



and 



3 



3Bin 7T 

7 ~ 2 



d0 c ou 



xl c , <9 C '» = 0, 



7T 



Furthermore, at 6>3 = 0, 

r 3 = x 2 Ax , xi = , x 2 = xf* 



and at 9 3 = ^ 



&Ay &By <5(7> 

r 3 = x 2 v , xi = r 3 y , x 2 — 
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We can now compute the derivatives of the 63 components of the full return map 
at and ^; we use the chain rule and make the assumption that the global parts of 
the maps only affect the derivatives by an 0(1) amount. We find that wc get different 
results, depending on the initial cross-section for the return map. This is consistent 
with the results derived in |17j using different methods. If we start on HJJ we have 

d6 A ™ _ f X?*, A '"=O 

Starting on H'J and HJ^ we have, respectively, 



■ 3 

„ V BY 



7T 



and 



7T 



2 



where 



VAX = & Ay — Sax-T^ H SaxSbx-, 

OBy ec 

VAY = 5 Ax — &Ay-T^- H 0~AyO~By, 

obx e c 

vbx = —r^- H $Bx + $a v SbxScx + (1 - Sx), 

"By e c 

8 X 

vby = -7-^ H S By + SaxSb v Scy + (1 - $y), 

obx ec 

vex = h Sa v Scx — SaxScx-t^-, 

ec o By 

A 8 By 
VGY = h dAxOCY — OAyOCY-^ • 

ec obx 
Note that the sign of the appropriate Vij determines the slope of the 63 part of the 
return map at #3 = or 63 = This in turn determines the stability properties of 
the invariant subspaccs at #3 = or 63 = -| in the full return map. 
The following relations hold between the constants defined above: 

vax8cx = v C x + —(8x - 1), v A y8cy = v C y H (8y - 1), 

ec ec 

Vbx8ax = VAX + (jr^- - lj (8 X - 1), V BY 8 Ay = V A Y + ~ y ~ !)' 



vcxSbx = vbx + 1 - t — (8x - 1), v C y8b v = vby + 1 - t — (8y - 1). 

V S ByJ V 8 Bx J 

If 8x is sufficiently close to 1, then vax, "m and vex all have the same sign; since 
we are interested in resonance phenomena for which Sx ~ 1, wc will assume this is 
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oA ia QA ia 

(a) VAX-.VAY > (b) u AX > 0, u AY < 

Fig. 4.1. Case I: 9% component of the full return map, starting on H^ 1 with a fixed value of 
X2, for two choices of signs of the i/ij . 



the case. Similarly, if Sy is sufficiently close to 1, then vay, vby and vqy all have the 
same sign; we will assume in the following that this is the case. This assumption means 
that the stabilities of the invariant subspaccs at #3 = and 63 = § arc independent 
of the section on which the composition of the return map starts. 

Away from resonance, it is possible that, for example, vax > and vbx < 0. It 
is precisely this type of condition which gives the very delicate stability properties of 
the subcycles of the network that is seen in [17]. There, a subcycle may appear to 
be attracting if nearby trajectories are observed as they pass through one Poincare 
section, but may seem to be repelling if trajectories are observed at a different Poincare 
section. This type of stability cannot be seen for objects such as periodic orbits or 
equilibria in flows. In this article, we only consider the case close enough to resonance 
when this phenomena does not occur, and hence need only consider composing the 
maps starting on HJ^ 1 . 

Consider first the case for which vax, vay > 0. Then, for fixed X2, the graph of 
9 Ain as a function of 9 Ai » has flat sections near 9 Ain = 0, f , 7r, that is, ^ A |° — » 
as X2 — > at these values of 9 Ain . Since almost all trajectories pass close to X, —X 1 
Y or — Y, almost all trajectories will return to H™ with a value of #3 approximately 
equal to 0, ^ it or As a consequence, the sections of the graph of 8 Aia between 
the flat sections will be steep. Figure 14.11( a) shows schematically the shape of the 
graph of 9 Aia as a function of 9 AiIL for trajectories with a fixed value of X2- 

The width of the small flat section near ^ can be computed. Points in this part of 
the graph correspond to orbits which pass close to Y , and the left boundary is given 
by the preimage of ^ — 9 tB under the map ^ ab where 9 tB was defined in (14. 4[) . 
We define 9 tA to be such that this preimage is ^ — 9 tA . 

We can compute 9 tA using the approximations of <f>A given in section 14.11 and 
assuming the global map has only an 0(1) effect. Since 9 tB <^ 1, we approximate the 

is 



4>A map as 



logr 3 = S A V log x 2 , 
tan(9 Bi » = tan6 3 x 6 2 Ay ~ SAx . 

Approximating tan# £B ~ l/9 eB (and similarly for 6 €A ), we have 

9 tA = 6 eB (x 2 ) Sav ~ 5a!1 ' 

6 By . 

— '3 x 2 




where we define 

<7 = O^y — OAs 7 • 

The right boundary of the fiat section near ^ can be found by symmetry, and so 

the width of the fiat section near ^ scales like x 2 Bt . 

Now we consider other cases of the signs of vax and vay ■ Note that 

vax = (J H ox and = — t — cr H ay. 

Thus, if a > 0, then i/^x > 0, while if a < 0, then i/^y > 0, and so it is not possible 
to have both vax < and vay < 0. This leaves the cases where vax and vay have 
opposite signs. If a > 0, then either v^y > 0, as considered above, or vay < 0, 
meaning that there is no small step near 63 = ? , and the invariant subspace #3 = ^ is 
repelling; a sketch of the component of the return map in this latter case is shown 
in figure l4~TT b). 

If a < 0, there are again two cases, similar to those described above, but with 
the roles of 9 = and 9 = -| reversed. We believe the dynamics for a < will be 
analogous to that for a > (only with this reversal) and so consider just the case 
a > for the remainder of this article. 

It is useful here to summarise the conditions we now have on the eigenvalue ratios 
in Case I. 

• By assumption A7, we have Sbx < fej- This implies that in a neighbourhood 
of Sx = Sy — 1 we must have 5c x 3 Ax > 5cY&Ay 

• We additionally choose to impose a > 0, and specifically want this to hold 
when Sx = Sy = 1- Since 

a > => SAySBy > SaxSbx 

we require Sex > Scy so that a > in a neighbourhood of Sx = Sy = 1- 
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• Together these conditions imply 



Sex 8 Ay 5bx 
Scy Sax Sb v 

Note that j^- could be greater or less than 1. 
Fixing 8bx and imposing the requirements < 1, > 1 and a > still gives 
us the freedom to vary both Sx and Sy above and below 1. 

4.2.2. X2 component. We now compose the three local and global maps start- 
ing on H™, using initial values of 9 AiTl and X2 with < #3 < J and X2 small, and 
focus on what happens to the X2 component of the map. We will eventually end up 
with an approximate map of the form 

X2 -> Dx 5 2 

where D and 8 are functions depending on 9 Aia and X2- In other words, the amount 
of contraction or expansion of the X2 coordinate of an orbit in one circuit of the 
network will depend on the initial condition for that orbit; this is a consequence of 
the network structure and is different to the case for maps modelling the dynamics 
near a single heteroclinic cycle. To capture this effect, in the following we write down 
the contraction or expansion rate of each the local maps as a function of the incoming 
coordinates for that local map, then rewrite the incoming coordinates as a function 
of the initial conditions of the orbit on H^ 1 . Thus, the functions we obtain for the 
contraction rates at B and C will depend on 9 Ai " and the value of X2 on H™. 

We use the approximate forms of the local maps derived in section I4z.ll making 
use of the assumed form of the dynamics at C. We will also assume that the global 
maps multiply the small variable by a ^-dependent order-one constant (as described 
in section . so r$ = DB{9s,)r A ° ut etc., and that the 63 parts of the global maps 
do nothing, that is, 9b(9^) = 6*3, and so 9 Bin = 9 AaMt etc.). This will give a distorted 
view of the correct picture, but the distortion will only be slight, since the dynamics 
is dominated by the local maps. 

We focus our discussion on the interval < 9 Ai " < ^ ; this can be extended to 2n 
by symmetry. To allow for this, we will include absolute values in expressions such as 
(for example) log | cos 8 AiIL \ . 

We divide the interval < 9 Aln < f into two regions, which are the different 

regions of validity of the approximate local maps 4>b and 4>c ■ The boundaries of the 

s b v 1 

regions depend on the value of r% on H'J as given in ()4.5|) . We have 9 CB = rg Bx 
and the two regions are given by cot 9 BiD 9 eB and cot 9 Bia <§C 9 eB . As computed 
in the previous section, the two regions can also be defined on H^ 1 as cot 9 Ain ^> 9 eA 

5 By 

and cot6> Ai " < 9 tA , where 9 tA = x^*" ■ 
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Region 1. First, consider the region < 9 Ain <C | - 9 tA , so ta,n9 Ai » x 



< 1. 



After 4>a and ^ab we have 

tan 9 Bi " 



tan 0^" x 



log r 3 =\ogD B + 5 a [0 A ' n ) log x 2 , 



where 



6 A (9 Ai ») = 8 Ax + 



log I cos 9 J 



I log ^1 + tan 
2 



2 nA iD Zi&Ay-SAx) 



log X2 2 log X 2 

Although 8a{6 Aiti ) depends also on log.T 2 , we omit specifying this dependence in the 
argument (and in i5s(6'" 4in ) and 5c (9 Aia ) below) to simplify the writing. 

Since we are in region 1, trajectories visit the X part of the second local map. 
Therefore, after <pB and $bc, we have: 

tan6» Cin =tan6» B » |r 3 cos s '" 
log Xi = log D c + S B (# Ain ) log r 3 



where 



s B (e J 



Sbx + $Bx 



log cos 9 



6 Bx 1 



log r 3 

^og^l + tan 2 0*** xl lBA "- SAm) 
2 5 A (9 A ^) log x 2 



Now, tan# in is small compared to 1, since we are in the region where cot# Sin 3> 



r 3 Bx . This follows from noting that 

, Sr. 



tan9 Cin =tan9 B '"r 3 By |cos0 B H 4 * 



tan 9 B "r* B * 



sin 9 



Bin 



3 * sgn (tan( 



The hrst term is small by assumption, and the second and third are at most 1 since 
5bx < 5 By (assumption A7). Therefore we use the X part of the map at C, and get, 
after <\>c and ^>ca- 

logx 2 = log D A + <5 c (6 |Ain )logxi 
where x 2 is the value of x 2 on H^ 1 after one full circuit of the network and 
Sc (e A *) = Sex - CC l X r CCY log (1 - tan 2 ^ 



= 5, 



cx 



2\\ogxi 

CCX — CCY 

2\6 A (9 A ^)S Bx \ogx 2 
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tan^6> 



2 ad, 



since tan 2 9 Cin <^ 1. 

Substituting for Xi in the above expression for X2, we find that 

logx 2 = logL> A - + S(9 A '") \ogx 2 , 

where 

\ogD x » Iog£U(0) + <5 C x (lo gj D c (0) + S Bx log£> B (0)) 



and 



( 1 log(l + tan 2 A -X2 (5 ' 4! '~ 5A * ) )\ 
X ~ 2 M^.»)loga- 2 J 

x ( S C x + ^tt t^t-^t; — ; tan 1 " 6»° ln 



2A* A (fr Ata )*Bxlogxa 



1 log (l + tan 2 A >* x f s ^- 5 ^ 



Sx + SbxScx- 



2 log x 2 

log J cos 9 Ain J 

fog CC2 



We have ignored the correction term in Sc(9 Ain ) since it is much smaller than that in 
Sa (0 a -)- 

In this region, S(9 Ain ) ranges between Sx (when 9 Ai " = 0) and S Ay SB y ScX: since 

>° 

at the edge of region 1, cos# ln ~ x 2 Bal ■ 

^JLa 

Region 2. Now consider the region with ta,n9 Ain x 2 Bx S> 1. Orbits with 9 Ain 
in this region will visit the Y parts of all three maps. Note that since 6 By > S Bx , the 
above assumption also implies that tan^" 4 '" x S 2 Ay <5 ' 4x 1. 

For <j>A in this region we write 

log |sin^»| 1 log (l + ^t 2 9 A * x 2 ^- 5 ^ 



s A {e A ^) = s Av + 



2 



log x 2 2 log x 2 



and we can approximate SA(9 Ain ) by 5 Ay- 
After 4>b, we find 



icot 2 B '"r? (1 ^ 

s B (e A ^) =5 Bv + s By ^- 



log r 3 



icot 2 9 A ^ 



Sa V fog X2 
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and after <fic we find 



6 c (e Ai ") = 5 CY 



2\8Ay$Byl0gX 2 



CCX — CCY 




) 



The corrections to the B and C parts of the map are small and comparable, 
but large compared to the correction to the A part of the map, so we find, for 

tan6» A '"x 2 Ba = > 1, 



In this region, the correction term could be of either sign since Sex > &cy- However, 



4.3. Case I: Resonance of a single subcycle. We can use the results derived 
in the previous section to consider resonance bifurcations of a distinguished subcycle 
within the Case I network. These results could be derived using the traditional cross- 
sections (as is done explicitly in [H]), and the results would be identical. However, 
rather than repeat that analysis, we show how these results can be achieved using 
our new methods. Specifically, we consider the subcycle of the network given by A — > 
B — > X — > A, which lies in the subspace j/3 = 0. This cycle cannot be asymptotically 
stable since B has a two-dimensional unstable manifold. 

The dynamics near this cycle are described by a two-dimensional map. Using the 
results of the previous section, it can be shown that the return map starting on H 1 ^ 
is given by 



If we start on a different section, the map will be similar, with, e.g., x 2 replaced by 
r3, and vax replaced by vby- 

The fixed point in this map at 63 = x 2 — corresponds to the heteroclinic cycle. 
We know the cycle cannot be asymptotically stable, but it can be attracting if Sx > 1 
and v>ax > (as discussed above). 

A resonance bifurcation of the heteroclinic cycle occurs when 6x = 1. This 
bifurcation creates a fixed point of the map at 6*3 = 0, x<x = D 1 ^ 1 Sx \ which is 
also in the subspace 2/3 = 0. Furthermore, it is straightforward to show that if 
< Dx < 1 then a periodic orbit occurs for Sx < 1 and so the bifurcation is 
supercritical, while if Dx > 1 then a periodic orbit occurs when 5x > 1 and so the 
bifurcation is subcritical. If this bifurcation occurs supercritically, then the resulting 
periodic orbit will be asymptotically stable. That is, we have the possibility that the 
resonance bifurcation is from a heteroclinic cycle that is not asymptotically stable 
but it produces a periodic orbit that is asymptotically stable. To the best of our 
knowledge, this scenario has not been reported before. 




in the limit of small x 2 , the value of S(8 AiI1 ) in all of region 2 is S Y - 



x 2 -> D x x 



Sx 



1 1 
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4.4. Composing the maps: Case II. Wc repeat the above calculations for 
the network with complex eigenvalues. There are again two regions, given by the 
same conditions as before. Due to the rotation of 63 at A, the regions are defined 
on Hg rather than H™, but we could map these back to H™ using the expression 
6 B " = 6 A *--^logx 2 . 

We again begin by considering the #3 components of the maps at # Bin = and 
6 Bin = t£. These points are not subspaces in this case (as they are in Case I), but 
can still give us information on the geometry of the #3 part of the return maps. The 
calculations proceed exactly as before, except that cax = ca v = ca- This means we 
have a simplification and find 



vax = Sa I 1 - 



JBy 



8a5l 



which must be positive. The relationships with vbx and vqx given in section 14.2.11 
imply that in addition vbx > and vqx > 0. Thus the 6 '3 part of the return map 
will have a small gradient close to the point where 9 Bia = 0. 

The computation of vay, vby and vcy follows as in Case I, and again we see 
that they all have the same sign so long as we are close enough to Sy = 1- Wc assume 
that this is the case, and further, that they are all positive, as before. The dynamics 
in the case where vay < is very similar. Thus, again, we need only consider the 
return map starting on H^. 

The graph of 9 Aia against 9 Ain will look very similar to that for Case I, shown in 
figure POT a) . except that as the initial value of x 2 varies, the graph will shift to the 
right or left. This is discussed in more detail in section [5721 

We next compute S(8 AiI1 ) for Case II, in exactly the same manner as for Case I. 
The only difference occurs after (f>A] now we have 5a{@ A>11 ) = 8a, which is a constant, 
and 



— log £2. 

&A 



The remainder of the calculations follow in exactly the same manner, and we find 
that in region 1, 



8{6 A ™) = 8 x + 8 B Jcx- 



and in region 2, 



log 


COS ^ 


6» A '" - — logx 2 ) 









log x 2 



S(0 a "') = Sy + 



•8cy 



(S C y - $cx)ec 
A 



cot 2 



J7 log x 2 ) X 



2 



2 log x 2 



5. Resonance of heteroclinic networks. We are now in a position to deter- 
mine the effect on the dynamics near each network of one or more of the cycles within 
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Fig. 5.1. The (Sx,&y) parameter plane, showing the definition of the quadrants A, B, C, and D. 

the network undergoing a resonance bifurcation. We focus on finding fixed points of 
the approximate return maps we have derived, which correspond to periodic orbits 
that make one circuit of the network before closing. 

Throughout this section, we start with a circle of initial conditions on H™ with 
fixed x-i and < #3 < 2w, and consider the values of xi and #3 when these trajectories 
first return to H™; we again refer to these values as £2 and #3, respectively. We use 
the approximations for the maps derived in section r4.ll to plot 'nullclines' of $3 and 
£2 on H™. A point (#3, x-i) € HJJ is said to be on the X2-nullcline (resp. ^-nullclinc) 
if the value of X2 (resp. #3) after one circuit around the network is unchanged (resp. 
unchanged modulo 27r). Fixed points of the Poincare map occur when the xi- and 
^-nullclines cross. We can identify these from the sketches of the nullclines, and are 
also able in some cases to identify the stability of the fixed points by considering how 
£2 and 63 vary close to the fixed points. 

We then discuss how the nullcline figures change as the quantities Sx and Sy are 
varied and are thus able to draw bifurcation diagrams. Figure [5~T1 shows the (Sx , Sy) 
parameter plane, and labels the four quadrants around the point Sx = Sy = \. In 
the following, we refer to these quadrants and also draw bifurcation diagrams as we 
traverse a small circle around the point Sx = Sy = 1. 

Recall that for both networks the return map has the general form 



where 0(63) is the constant arising from the global parts of the map and S (63) (which 
depends on X2 as well as #3) was calculated in section |U If 5(83) > 1 for all 6*3 and 
D(6j,) < 1 for all 63, then log £2 < log £2 for all #3 and all small X2- Hence, the 
network is asymptotically stable. If (5(6*3) > 1 but 0(63) > 1 for some 6*3, then for 
sufficiently small X2 , log £2 < log X2 and the network is still asymptotically stable. 
However, if 0(63) > 1 and X2 is large enough that log £2 > log-D^)/^ — S(9 3 )), 
then log £2 > log £2 and trajectories move away from the network. Thus, in the case 
that 0(83) > 1 for some #3, the basin of attraction of the network could be quite 
small as 3(63) — > 1 from above. 



logx 2 = log£>(6> 3 ) + S(9 3 ) log x 2 



(5.1) 
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Fig. 5.2. The figure shows 83-nullclines on H^ 1 (in blue) for the network in Case I, when 
vax > 0, and 8(83) > 1, so the network is attracting (guadrant A in fiaure \5.lV . The orange arrows 
indicate that X2 is decreasing, and the blue arrows indicate the direction of change of 83 . 



For simplicity, we thus consider only the case when D{6^) < 1 for all #3. This 
means the network is attracting and has a large basin of attraction if 6(03) > 1 for all 
#3, which makes it simpler to study what happens when 6(83) goes through 1. This 
condition on 0(83) is similar to assuming a supercritical bifurcation in other types of 
bifurcation. 



5.1. Case I: Computing nullclines. We first consider computing the #3 and 
X2 nullclines for Case I, the network with real eigenvalues. 



5.1.1. #3-nullclines. We begin by finding fixed points of the #3 part of the map, 
and using this information to draw #3 -nullclines on H™. Figure |4~T1 shows the value 
of 83 after one excursion around the network. There are fixed points at 83 = 0, w, tt, 

If vay > 0, then there are four further fixed points either side of \ and These 
additional points are at (approximately) -| ± 8 eA and ^ ± 8 eA , and so get closer to 
§ and ^ as 22 decreases. 

Figure l5?2l shows a sketch of the ^-nullclines in the case vay > 0. The distance 

from the curved 8- nullclines to -| scales like x 2 Bx ■ The blue arrows in the figure 
indicate how 83 changes under iteration of the map. This shows that the nullclines at 
83 = 0, §, 7r, ^ are attracting, but the additional (curved) nullclines are repelling. 
In the case that Uay < 0, the additional nullclines are not present and the nullclines 
at 77 and are repelling. 

5.1.2. S2-nullclines. We next construct the 22-nullclines. Our calculations are 
done explicitly for the region < #3 < ? but results for the remaining values of #3 
follow from symmetry. 
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Fig. 5.3. Case I: a sketch of 6(63) against 63 for fixed X2- 



The return map has the form given in (|5.1j) with 



5x + 5 Bx 5cx 



log I cos 3 1 



O<0 3 < f -0 £A , 



log X 2 



s(0 3 ) = { 



(8 C Y - S C x)ec\ cot 2 9 3 x 2 




Sy + 5 By ScY + 



A / 2 log x 2 



\-9 f - A < 3 < |. 



Figure IOI shows a sketch of 5(03) against 03. As discussed in section |4~T| in region 1, 
5(03) varies between Sx and 5m, where 5m = 5a v 5b v 5cx ■ Note that 5m is greater 
than both 6x (since a > implies that 5A y 5By > 5a x 5bx) and 5y (since 5cx > 5cy)- 
The existence of this maximum of 5{9 3 ) close to 5m is persistent in the limit of small 



could be local minima or maxima, depending on the sign of the factor in front of the 
second term in 5{9 3 ) in region 2. However, this correction term is much smaller than 
the correction term to 6(63) in region 1, and the value of 6(63) on the boundary of 
region 2 tends to 5y in the limit of small 

Thus, the maximum value of 5{9 3 ) is 5m, and the minimum, in the limit of small 
X2, is either 5x or 5y- In [TS], we showed that if min^ 5(9^) > 1, then the heteroclinic 
network is asymptotically stable, and if maxe 3 <5(0 3 ) < 1, then the heteroclinic network 
is completely unstable in that the basin of attraction has measure zero. Therefore, we 
expect to see stability changes, or resonances, of the heteroclinic network when 5x, 
5y or 5m pass through 1. 

Intuitively, we expect to find fixed points near 03 = if Sx < 1 (but close to one) 
and fixed points near #3 = ^ if <5y < 1 (but close to one). To check this, we find the 
;E2-nullclines explicitly by finding solutions to the equation 



x-2- Note also that 03 = 0, it are always local minima of 5(9%) but 0; 




3tt 
2 



logz 2 = logD(0 3 ) + 5(9 3 ) log x 2 . 
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If such solutions exist in the region 8 < % — 8 eA , then we have 



logx 2 = log -D(0 3 ) + S x \ogx 2 +5 Bx 8 C x log | cos6» 3 | 
which, after rearranging gives the curve describing the nullclines: 

logx 2 = - — (log-D(0 3 ) + £ B xfcxlog| cos0 3 |) • 

1 — Ox 

Since we assume D{9^) < 1. we require 5x < 1 for solutions in this region, as expected. 
This curve has a maximum at 83 = 0, where log £2 = log.Dx/(l — 5x)- For later 
convenience, we define x x = 

Suppose now that solutions exist in the region ~ — 8 eA < 63 < ^ . These solutions 
satisfy 

(8 C Y -S C x)e c \ cot 2 9 3 x 2 Bi 



log£ 2 = log D(8 3 ) + Sy log-T 2 + ySsyScY 
To leading order, we can write this as 



A / 2 log x 2 



logx 2 = - 1 log-Py, 

1 — Oy 

and hence for solutions in this region we require 5y < 1, as expected. For later use, 
we define x Y = Dy'" 1 ~ 5Y \ 

If Sy < 1 and 6m > 1, then there will be additional solutions at the boundary of 
the two regions, that is, where 83 ~ ^ — 8 £A , for x 2 < x Y . Note that the a;2-nullclines 
concerned have the same scaling (in terms of distance from ^ ) as the additional $3- 
nullclines (which exist only if v&y > 0). Thus, to determine the relative positions of 
the two sets of nullclines, and to work out where the nullclines cross, we would have 
to include more details about the global constants. In practice, it is likely that both 
cases are possible; we discuss the possibilities further below. 

In figures 15.41 15.51 and 15.61 we show sketches of the #3- and X2-nullclines in the 
quadrants B, D and C around the point 5x = Sy = 1, sufficiently close to that point 
so that 8m > 1- We show figures only for the case vay > 0, and so the additional 
£2-nullclines are present, but discuss the case vay < below. 

In figure HT4I Sx < 1 and 6y > 1, and we can see that a stable fixed point occurs 
at 83 = 0, x 2 = x* x (and similarly at #3 = 7r, by symmetry). In figure f5T6l 5y < 1 
and 5 X > I. In the case shown, vay > 0, and there is the possibility of cither one 
or three fixed points appearing close to #3 = | at resonance (and also near 4^, by 
symmetry). The figure shows the case where the additional X2-nullclines lie further 
from ^ than the additional #3-nullclines, and three fixed points are created, one stable 
and two of saddle type. A second possibility is that the X2-nullclines lie closer to \ 
than the #3-nullclincs and there is only a single stable fixed point created as by passes 
through 1. If vay < 0, then there would also only be a single fixed point created as 
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Fig. 5.4. The figure shows nullclines on H'jj 1 for 83 (blue) and X2 (red) for Case I, in quadrant 
B of the (8x , 5Y)-plane. The orange and blue arrows denote the direction of change of X2 and 83 
respectively. The pink dots indicate stable fixed points of the map. 



5y decreases through 1, but in this case it would be of saddle type as the nullcline at 
$3 = § would be repelling. 

In figure 1531 5x,5y < 1, and we show the figure for 5 m > 1- Both sets of fixed 
points described above exist, and all the nullclines continue to exist as loga^ decreases 
to —00. In this case, the fixed points created in the two resonance bifurcations at 
5x = 1 or 5y = 1 do not interact with each other, a consequence of the two red 
nullclines being distinct from one another for arbitrarily small x^- 

Finally, in figure [57H we show the case where 5 m < 1. This is still in quadrant C of 
the (Sx, Sy) plane, since 5m > Sx,5y. In this case the a^-nullclines created in the two 
resonance bifurcations of the individual sub-cycles have joined up, and rz^-nullclincs 
exist only for a finite range of log xi . The additional resonance bifurcation that occurs 
when 5m passes through 1 has the possibility of creating further fixed points near the 
additional ^-nullclines, if they exist (i.e., if vay > 0), and if they were not already 
created in the 5y = 1 resonance. 

5.1.3. Bifurcation diagrams. We now use the nullcline sketches to draw bi- 
furcation diagrams. In figure [5T8l we show a bifurcation diagram obtained as a circle 
is traversed around the point 5x = 5y = 1 in the (5x, Sy) plane. We assume we are 
close enough to this point so that 5m > 1 and hence that the periodic orbits created 
when 5x passes through one are not connected to those that arise when 5y passes 
through one. 

If vay > 0, there are two cases to consider: cither the a^-nullclincs arc closer to 
^ or further from ^ than the curved #3-nullclines. In the first case, the only equilibria 
are at 03 = and 9s = § • In the second case, there will be further equilibria; one 
possibility is shown in the right panel of figure 15.81 

The supplementary online material contains a movie showing how the nullclines in 
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Fig. 5.5. The figure shows nullclines on H'J 1 for 63 (blue) and X2 (red) for Case I, in quadrant 
C of the (5jf,<5y) plane. All lines, curves and dots have the same interpretation as in figures \ 5.4\ 
except that the green dots indicate saddle fixed points. 



loga; 2 
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Fig. 5.6. The figure shows nullclines on H'jj 1 for 63 (blue) and xi (red) for Case I, in quadrant 
D of the (Sx,&y) plane, in the case 6~m > 1 All lines, curves and dots have the same interpretation 
as in figures \5.4\ and 15. 51 



Case I vary as a circle of radius 0.02 around 5x = &y = 1 is traversed in the (Sx,5y) 
plane. In this movie, we kept cr fixed at 0.05, = 0.07, = 0.93, fz* = 1.08, 
and chose the other coefficients such that the value of for which changes 
from a local minimum to a local maximum is 0.07. The red solid curves are the 
small nullclines in region 1, the green solid curves are the small nullclines in region 2, 
and the blue solid curves are the 6*3 nullclines. The red and green dashed curves are 
the approximate small nullclines computed above. Regions 1 and 2 are separated by 
green dashed curves. As the point Sx = Sy = 1 is circled, the region 1 and region 2 
small nullclines appear and disappear as the lines 5x = 1 and 8y = 1 are crossed 
respectively, leading to the creation or destruction of fixed points near X or Y. 
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Fig. 5.8. The figures show the creation of periodic orbits as a circle is traversed clockwise 
in (5x,$y) space around the point &x = Sy = lj for the Case I network. The left hand figure 
shows the case when the x^-nullclines are inside the curved d^-nullclines, and the right hand figure 
shows the case when the X2-nullclines are outside these 63-nullclines. The labels A, B, C and D 
correspond to the quadrants labelled in figure I5.il The pink curves indicate stable periodic orbits 
and the green curves indicate saddle periodic orbits. (Note that the pink curve actually represents 
two symmetry-related orbits and the green curve four). 



5.2. Case II: constructing nullclines. A similar analysis can be performed 
for the network with complex eigenvalues. 

5.2.1. #3-nullclines. We will assume that Vay > 0; the situation for vay < 
in Case II has only very minor differences. 

Plotting the value of 63 as a function of 6*3 for some fixed initial value of X2 gives 
a schematic picture similar to that shown in figure l4~IT a). However, differences are 
noticed as the value of X2 is decreased. Specifically, the effects of reducing the initial 
value of X2 include those given above for Case I, i.e., the steep portions of the graphs 
become steeper, and the small 'step' becomes smaller, but the additional time spent 
in a neighbourhood of A when X2 is smaller means that the value of #3 is 'rotated' for 
longer due to the complex eigenvalues (specifically, 6* Sin = 9 Ai " — ^- log £2). This has 
the effect of shifting the graph of #3 to the left as X2 is decreased. This means that 
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Fig. 5.9. Case II: nullclines of 63 (blue curves). Compare with figure \4- l\ which shows how 83 
varies at fixed X2- The blue arrows indicate the direction of change of 6. 



the topology of the #3-nullclincs is different in Cases I and II, as we now explain. 

For the value of 22 shown in figure |4~T1 there are four points at which the value 
of #3 is the same after one circuit of the network. These points are thus on the #3- 
nullclines. As x 2 decreases, the graph of 9 3 moves to the left, and thus the four 'fixed 
points' in the #3 map come together and disappear in pairs, in a manner similar to a 
saddle-node bifurcation in a map. There are then some values of xi for which there 
are no fixed points in the #3 map. If X2 decreases so that the value of #3 — ^ logX2 
has changed by 2tt, then the graph in l4.1l will have rotated back to its original position 
(except that since x-i will now have decreased, the vertical parts will be steeper and 
the small step smaller, as discussed previously). 

Figure [53] shows the location of the £>3-nullclines on H 1 ^ . The vertical gap between 
the nullclines is such that the difference in log 22 is 27r ^ A ■ Note that is a cylinder, 
and each of the # 3 -nullclines is topologically a circle around the cyclinder. There is 
an infinite number of these nullclines. The larger approximately vertical portions of 
each #3-nullcline should appear at #3 = and #3 = n, by our assumption that the 
global parts of the #3 maps do nothing. However, for clarity, in this and following 
figures we show these portions of the curves slightly away from and n. This has no 
effect on the topology of the intersections of the ^-nullclines with the cc2-nullclines 
we describe later. 

In figure |5"U1 we also show how #3 changes away from the nullclines, marked with 
blue arrows. We only include these close to the nullclines, as since #3 is a circular 
variable, it does not make sense to say whether 63 is increasing or decreasing when it 
is changing by a large amount. Thus the direction of change of #3 can change from 
right to left without crossing a nullcline. 

5.2.2. a^-nullclines. Determination of the existence and shape of the a^-null- 
clines proceeds exactly as in Case I, except for consideration of additional rotation as 
X2 decreases, as for the #3-nullclincs. Thus, the a^-nullclines for Case II will look the 

32 




same as in Case I except that the 83 coordinate is replaced by #3 — ^- log X2- In other 
words, the $3 coordinate of the nullclines rotates to the left as X2 decreases. 

Figures f5. 101 15.11] and !5.12l show the 83 and a; 2 -nullclines for quadrants B, C and 
D of the (5x,5y) plane respectively, for 5m > 1. In these cases, the X2-nullclines 
exist for arbitrarily small X2, and so there will be an infinite number of intersections 
of the 83- and X2-nullclines, and hence an infinite number of fixed points in the map 
or periodic orbits in the original flow. 

In figure I5.10| 6x < 1 and <$y > 1. As 6x decreases through 1, fixed points 
are created in saddle-node bifurcations for 83 m 0, n and with X2 ~ x x . In each 
saddle-node pair, the larger amplitude solution is initially stable, and the smaller is 
of saddle-type. As 5x changes, it is likely that these fixed points undergo period- 
doubling or other types of bifurcation, and hence their stabilities may change. 

In figure I57H21 6y < 1 and 5x > 1- As 6y decreases through 1, fixed points are 
now created in saddle-node pairs near 83 « f , ^ and with X2 ~ x\. Again these 
fixed points will initially be created in stable-saddle pairs, but due to the small step 
in the 83 map and the shape of the # 3 -nullcline, we expect the 83 coordinate of these 
points to change rapidly as 5y is varied, and expect some of them to undergo stability 
changes too. 

Figure IBTTTl shows the situation when 5x,5y < 1, 6m > 1; as in Case I, sets of 
periodic orbits from the resonances at 5x = 1 and 5y = 1 co-exist in this quadrant. 
Finally, in figure [5TT31 we show the case 5m < 1- Here the a;2-nullclines only exist for 
a finite region of loga;2, and hence there are only finitely many fixed points. Thus, 
the resonance bifurcation which occurs at 5m = 1 in the complex case results in the 
disappearance of infinitely many periodic orbits. 
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5.2.3. Bifurcation diagrams. Figure 15.141 is a bifurcation diagram showing 
how periodic orbits are created and destroyed as a circle is traversed around the point 
fix = fiy = 1) assuming that 5m > 1- 

The supplementary online material contains a movie showing how the nullclines in 
Case II vary as a circle of radius 0.02 around Sx = fiy = 1 is traversed in the (fix, fiy) 
plane. In this movie, we kept a fixed at 0.05, = 0.07, = 0.93, = 0.5, 
and chose the other coefficients such that the value of for which S(^) changes 
from a local minimum to a local maximum is 0.07. The red solid curves are the small 
nullclines in region 1, the green solid curves are the small nullclines in region 2, and the 
blue solid curves are the 83 nullclines. Regions 1 and 2 are separated by green dashed 
curves. As the point fix — fiy = 1 is circled, the region 1 and region 2 small nullclines 
appear and disappear as the lines fix = 1 and fiy = 1 arc crossed respectively, leading 
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Fig. 5.13. Case II with 8y < 1> Sx < 1 (quadrant C), and with Sm < 1: nullclines for 9s 
(blue) and x 2 (red). Dots and arrows have the same meaning as in figure [5. \ 0\ 
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Fig. 5.14. Case II: bifurcation diagram, showing the creation of periodic orbits as a circle is 
traversed clockwise around the point Sx = Sy = 1 in (Sx,Sy) space, where Sm > 1 on the entire 
circle. The labels A, B, C and D correspond to the quadrants shown in figure 15. 1\ The pink and 
green curves correspond to the fixed points coloured pink and green in figures \5.1 Olf/rTgl stabilities 
may change along these curves. 



to the creation or destruction of infinite numbers of fixed points. 

5.2.4. Chaotic attractor. It was noted in [TH] that chaotic attractors can be 
found close to the Case II network when Sx < 1 and 6y > 1; it was argued that 
trajectories passing near X would be pushed away from the network (since 8x < 1) 
while trajectories passing near Y would be pulled towards the network (since Sy > 1)- 
A balance between contraction and expansion for orbits that pass repeatedly near X 
and Y could then be achieved, and may result in chaotic dynamics. 

Here we refine this argument, supposing first that we have a chaotic attractor, 
and then looking more carefully at the conditions needed to allow it to exist. This 
hypothesized chaotic attractor will have a range of values of logX2 on H™, and so 
there will be a corresponding range of values of 9 AaMt = 6 Airi — log Xi . If the 
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chaotic attractor is close to the network, then the range of 9 Aout will exceed 2ir, and 
could be many times 2n. In this case, orbits on the attractor will experience an 
overall contraction (towards the network) that is the average of S(0 Airi ), as given in 
section l4~4l We can approximate the average as: 



— / d(0)d0^dx-\ / log cost/ atH . 

71" Jo log X 2 7T Jo log X 2 



Note the inclusion of log D (the average of the global constant) in this expression. The 

contribution from the Y part of the cycle will be proportional to x 2 Bm , which is small 
compared to the 1/ log 22 term, and so has been dropped. The integral evaluates to 
— § log 2, so we find 



j r . logD - 5 Bx Scxlog2 

~ Ox H : ■ 

log x 2 

Finding x 2 so that 6=1 gives the expected distance of the chaotic attractor from the 
network: 



log D - S B Jcx log 2 
\ogx 2 ^ l~5x ^ ^ 

suggesting that the chaotic attractor bifurcates from the network at Sx = 1 in the 
same way as the periodic orbits shown in figure [5TT0l The term — Sb x 5cx log 2 is 
negative, which suggests that the chaotic attractor will be closer to the network than 
the periodic orbits. This issue is explored numerically in more detail below. 

Replacing the actual trajectory by the average in this way implicitly assumes 
that the distribtution of Amt is uniform. This will be a better approximation if the 
chaotic attractor is closer to the network, or if lu is larger. However, a non-uniform 
distribution would just lead to replacing log 2 by a different order-one number. 

Note that this estimate for the location of the chaotic attractor created in the 
5x = 1 resonance is independent of 5y , in contrast to the explanation offered in |18j . 

5.2.5. Numerical example. In this section we give an ODE that has a network 
of the type we are considering in this article. We give an example close to Sx = Sy = 1 
where there are a large number of stable periodic orbits coexisting with three chaotic 
attractors at the same parameter values. The equations are similar to those presented 
in pi]: 

x\ = X\(l — x\ — 2^2), 

x 2 = x 2 {l - x 2 2 - (1 + 5cx)xl - (1 + 6cy)Vs)> 

±3 = x 3 (l - (1 + 5a)x\ + s ^ Bx ^j x i - x i - C 1 + - uy-ix\, 

lis = 2/3 - (1 + 5 A )xj + (^^) + \)xl - ylj + uxzx\. 
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Fig. 5.15. Intersection of 16 different trajectories with the Poincare section Hj (h = O.OlJ, 
with Sx = $Y = 0.99. Reading from the top down, there are 11 stable periodic orbits, each separated 
in logarithm by 27r J A . Next, there are two period- doubled periodic orbits, then (between the horizontal 
lines) there are two distinct regions of chaos that visit only +X . Below the third horizontal line 
there is a third region of chaos that extends to —X as well; the average value of log r^ on this chaotic 
attractor is —166 (the average value of log 2:2 on H'^ is — 230). Black asterixes (blue plusses) indicate 
that the trajectory visits X (—X ) immediately after leaving the Poincare section. The boundaries 

of the cuspoidal regions are tan #3 = r 3 Bx ; inside these cusps, trajectories would visit Y or —Y. 



These ODEs have the fixed points A at (1,0,0,0), B at (0,1,0,0), X at (0,0,1,0) 
and Y at (0, 0, 0, 1). The constants 5a, 5b x , etc. are eigenvalue ratios with the same 
meaning as used throughout this article. 

We have carried out computations in each of the four quadrants indicated in 
figure 15.11 but present only one example here, for 5x = 5y = 0.99 (quadrant C in 
the (5x,5y) plane). The other parameters are 5a = 0.7143, 5b x = 1.4, 5 By = 1-5054, 
5 C x = 0.99, 5 C Y = 0.9207, A = 0.07 and ui = 0.5. The combination 5 M is 1.0645, a 
is 0.05, and all the u's are positive. The numerical methods are as described in [15] . 

In this example, the network is unstable and trajectories that start very close 
to the network move away from it. We have found 11 stable periodic orbits in the 
locations that would be expected from the considerations in section [5.2.21 Closer to 
the network than these, there are two period-doubled orbits and three distinct regions 
of chaos. The closest of these to the network has a reasonably uniform distribution 
of 6 in , and equation (|5.2j) is satisfied if we take the value of logD to be —1.34. The 
other two chaotic attractors have non- uniform distributions of Bin . We would expect 
there to be (possibly stable) periodic orbits that visit Y (since Sy < 1), but we have 
been unable to find these. Even if the orbits were stable, we would expect them to 
have small basins of attraction. 

The behaviour observed for parameters in quadrant B of figure 1531 (for example, 

37 



5x = 0.99, 5y = 1.01) is the same as that seen for quadrant C; since we were unable 
to locate periodic orbits that visit Y in quadrant C we do not notice their (predicted) 
absence in quadrant B. The behaviour in quadrants A and D (for example, Sx = 1-01 
and 5y = 1-01 or 0.99) is as expected from [18]: the network is attracting, and 
trajectories that start close enough to the network go towards it, repeatedly and 
irregularly switching between +X and —X, even though in region D, the network is 
not asymptotically stable (since dy < 1). In both regions, there are stable periodic 
orbits further away from the network. 

6. Resonance bifurcation of a single heteroclinic cycle with complex 
eigenvalues. To put in context the results we have found for resonance of our Case II 
network, it is helpful to look at resonance of an isolated cycle in which the linearisation 
of the vector field has a pair of complex conjugate eigenvalues at one equilibrium of 
the cycle. The cycle we consider is the same as one of the subcycles of the Case II 
network with itinerary A — > B — > X — > A except that at equilibrium B there is only 
a single positive eigenvalue, and hence, the unstable manifolds of all the equilibria in 
the cycle are one dimensional. Since we are interested in this section in orbits that 
lie near a single heteroclinic cycle rather than in a continuum of heteroclinic cycles, 
we can use much simpler forms for the local and global maps than in our analysis 
of the Case II network, and we are able to compute the full return map with ease; 
our analysis is analogous to that used for investigation of homoclinic bifurcations of a 
saddle- focus in, for instance, [T4l[T5] . Furthermore, existence and stability of periodic 
orbits near the cycle can be deduced from analysis of a single return map; there is 
no need to look at return maps defined on cross-sections near all the equilibria. We 
find that at resonance of this cycle an infinite number of periodic orbits appear in 
saddlc-node bifurcations, in a similar way to that seen for resonance in the Case II 
network. 

Specifically, we consider a system of ODEs in K 4 that is equivariant under the 
symmetries Ki, K2 and K3 as defined in (|2.1[) . (|2.2j) and (|2.5j) . and suppose that there 
are equilibria, £1, £2 and £3 on the positive x\, X2 and x 3 coordinate axes, respectively. 
These play the role of A, B and X. We assume that there is a connection from £x to 
£2 in the invariant (x±,X2) plane, a (single) connection from £2 to £3 in the invariant 
subspace defined by x\ =0 (this connection is not assumed to lie in a coordinate 
plane) and a connection from £3 to £x in the subspace defined by X2 = 0. The 
existence of invariant hyperplancs allows us to consider just the region of phase space 
where x\ > and X2 > 0. To simplify the discussion, we will also consider only 
trajectories that leave £2 with x 3 > 0, that is, we do not consider trajectories that 
visit —£3. 

The flow linearised about £x is given by 

mi = -riiti, x 2 = ei.T 2 , ±3 = -C1X3 - u>y 3 , y 3 = ujx 3 - ciy 3 , 
where fx, ex, c\ and uj arc positive constants, and where the u\ coordinate is obtained 
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from x\ after translation to move £1 to the origin of the local coordinate system. 
Near £ 2 , we use local coordinates u 3 and v 3 that are linear combinations of the global 
coordinates x 3 and y 3 , and local coordinate u 2 that is a translation of x 2 . The 
coordinate x\ is the usual global coordinate. The flow linearised about £2 is then 
given by 

X\ = -c 2 Xi, u 2 = -r 2 u 2l U3 = e 2 u 3 , «3 = ~s 2 v 3l 

where r 2 , e 2 , c 2 , s 2 are positive constants. The flow linearised around £ 3 is similar: 

x\ = e 3 xi, x 2 = -c 3 x 2 , u 3 = -r 3 u 3l y 3 = -s 3 y 3 . 

Here we use a translated u 3 coordinate but the other coordinates are just the global 
coordinates. 

It is convenient to use planar cross-sections near each equilibrium, For instance, 
we define 

H i n = {(ui,x 2 ,x 3 ,y 3 ) I \m\ <h, x 2 = h, \x 3 \ < h, \y 3 \ < h} 

and define H™, H 2 n , H3 11 and H 3 n in a similar and obvious way. We define cross-section 
H\ n slightly differently: 

Hf = {{ Ul ,x 2 ,x 3 ,y 3 ) I |«i| <h, 0<x 2 <h, xoe"^^ < x 3 < x e^ /u} , y 3 = 0} 

where the positive constant xq is chosen so that the heteroclinic connection from £3 
to £1 crosses at x 3 = xq, and the bounds on x 3 ensure that there is just a single 
intersection of the connection with the cross-section. 

Using these coordinates and cross-sections, it is straightforward to derive local 
and global maps. To lowest order, these are: 

01 (ui, x 2 , x 3 , 0) = ^ui (^j ei , h, x 3 ei cos log , 

(X 2 \ 7J . / Ul (X 2 

<p 2 (h, u 2 , u 3 , v 3 ) = (^j " 2 , u 2 (j^J " 2 , h, v 3 

(j) 3 {xi, h, u 3 , y 3 ) = (h, h (^j " 3 , u 3 " 3 , y 3 [^j- 

*i 2 (iii, h, x 3 , y 3 ) = (h, e 2 , ax 3 + by 3 , cx 3 +dy 3 ), 
tyjsfcl, u 2 , h, v 3 ) = {fxi, h, e 3 , gxf + ju 2 + kv 3 ), 
^3i{h, x 2l u 3 , 2/3) = (ei, mi 2 , x + nu 3 + py 3 + qxj, 0), 

where e,, a, b, c, d, /, g, j, k, m, n, p and q are constants. Composing these maps in 
order gives the return map / : H" 1 — > Hj n , which to lowest order is: 

x 2 , x 3 , 0) = fei, A x x 2 ^e 3 cos ^A 2 - ^- logx 2 ^ , A 3 , oj , (6.1) 
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where A\, A 2 and A 3 are constants and 6 = {cic 2 c^) / ■ This map is defined 
for sufficiently small \u\\, x 2 and X3, with x 2 > and x 3 > 0. In addition, the map is 
only defined for values of x 2 for which the cosine is positive. 

At lowest order, fixed points of the return map occur for m = ei, £3 = A3 and 



where A = AiA^ 2 " 3 > 0. Equation (|6.2[) is very similar to the type of fixed point equa- 
tion obtained in analysis of a Shil'nikov homoclinic bifurcation in a non-symmetric 
context |T4j[15] , with the differences being that (|6.2[) has an exponent on the cosine 
term and no bifurcation parameter on the left hand side of the equation; this last 
difference reflects the fact that we are interested in bifurcations that occur as 6 varies 
and the cycle persists but passes through resonance rather than as the cycle is created 
or destroyed by relative movement of its stable and unstable manifolds. 

Figure IBTT1 shows schematically graphs of the functions h\(x 2 ) = x 2 and 



for qualitatively different choices of 6 and A; fixed points of I correspond to inter- 
sections of these two graphs. As can be seen in panel (a), if 6 < 1 there will exist 
infinitely many fixed points of the return map, with the fixed points accumulating 
on the origin. This corresponds to the existence of infinitely many periodic orbits 
accumulating on the heteroclinic cycle. On the other hand, as shown in panel (b), if 
6 > 1, there will be no fixed points of the return map in the vicinity of the origin; this 
corresponds to there being no periodic orbits lying in a sufficiently small neighbour- 
hood of the heteroclinic cycle. The situation for the case 6 = 1 depends on the size 
of A; if A > 1 we expect infinitely many periodic orbits to exist when 5 = 1, while 
if A < 1 there will be no periodic orbits in a sufficiently small neighbourhood of the 
heteroclinic cycle when 6 = 1. 

Consideration of the possible transitions between the different cases shown in 
figure 16.11 now enables us to sketch schematic bifurcation diagrams showing the be- 
haviour of periodic orbits near the resonance bifurcation. As shown in figure 16.21 in 
the case that A > 1, for sufficiently large 5 > 1 there will be no periodic orbits in a 
small neighbourhood of the heteroclinic cycle. As 6 decreases, periodic orbits will be 
created in pairs in saddle-node bifurcations, with the saddle-node bifurcations accu- 
mulating on 6 < 1 from above, thus producing an infinite number of periodic orbits 
for all positive 6 < 1. For the case A < 1, there will similarly be no periodic orbits 
near the heteroclinic cycle for sufficiently large 6 > 1 and infinitely many periodic 
orbits for 6 < 1, but the periodic orbits now appear on the opposite side of the res- 
onance bifurcation; an infinite number of saddle-node bifurcations of periodic orbits 
accumulate on 5 = 1 from below, so an infinite number of periodic orbits will appear 
all at once as 6 decreases through 1. 




(6.2) 
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Fig. 6.1. Schematic diagrams showing the location of fixed points of the return map I, equation 

fO) . for different choices of S and A = AiA^ 2 ' 3 . (a) 5 < 1; (b) 8 > 1; (c) S = 1, A > 1: 
(d) 5 = 1, A < 1. Each panel shows the relative position of the graphs of h\(x2) = X2 (in red) and 

^2(2:2) = Ax^ {cos \ A.2 — ^-logX2^ " 2 " 3 (in blue). The black dotted curves show the graphs of 
hi{x2) = ±Ax^- Fixed points of I correspond to intersections of the graphs of hi and /12. Note that 
the shape of the graph of /12 where it cuts the X2 axis will depend on the exact value of the exponent 
°3°3 . For the purpose of illustration, we have drawn the case fjftj- = 1, ctnd included hi(x2) below 
the X2 axis to make the graph easier to read. We are not concerned with values of X2 for which /12 
is non-positive or undefined. 



Approximate 5 values for which saddle- node bifurcations of periodic orbits occur 
can be computed by comparing the graphs of hi(x2) and /12 plotted in figure UTTl 
Specifically, making the approximation that saddle-node bifurcations occur at X2 val- 
ues for which h 2 has a local maximum allows us to compute that, to first order, 
successive saddle-node bifurcations occur at 

5 n = l< ul ° gA 



ei(2nn - A) ' 

from which it follows that the saddle-node bifurcations accumulate on d = 1 exactly 
as derived schematically in the previous paragraph. We have not computed the values 
of S for which the node-type periodic orbits created in each saddle-node bifurcation 
are stable, but note that these nodes will likely change stability in period doubling 
bifurcations near the saddle-node bifurcations, and may undergo cascades of period 
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(a) 



(b) 



Fig. 6.2. Schematic bifurcation diagrams for the example of an isolated heteroclinic cycle with 
complex eigenvalues, showing periodic orbits that occur near the resonance bifurcation at 8 = 1. (a) 
A < 1; (b) A > 1. Stability of the periodic orbits is not indicated. 

doubling bifurcations leading to chaos, just as occurs in homoclinic bifurcations of 
saddle- foci [T4j[15] , and indeed as suggested by the numerical results in section 15.2.51 
The bifurcation diagrams obtained for resonance of this single cycle are completely 
consistent with the bifurcation diagram for resonance of our Case II heteroclinic net- 
work; compare hgures [5".14l and [p r ?2T a). This leads us to conjecture that the appearance 
of infinitely many periodic orbits near resonance of the Case II cycle is primarily due 
to the complex eigenvalues in the network not to the network structure. We note, 
however, two points. First, the analysis in this section explicitly requires that all the 
equilibria in the network have one-dimensional unstable manifolds and so, while our 
results are suggestive, they do not apply directly to the network example. Second, 
our analysis of the Case II network focusscd on periodic orbits that made just one cir- 
cuit of the network before closing and therefore excluded orbits that explored much 
of the network structure. It is likely that the bifurcation diagram for the network 
example contains sequences of saddle-node bifurcations additional to those we found. 
For instance, there might be infinite sequences of bifurcations producing orbits that 
make one or more visits to Y interspersed with visits to X. Such bifurcations could 
be regarded as arising from the network structure; investigation of this possibility is 
left to future work. 

7. Discussion. In this article, we have investigated resonance bifurcations in 
two robust heteroclinic networks; we believe this is the first time any examples of 
network resonance have been systematically studied. The networks of interest have 
both previously been studied [T71IIH], and consist of a finite number of equilibria 
connected by heteroclinic connections. An important feature of both networks is 
that several of the equilibria have two-dimensional unstable manifolds, which results 
in the existence of an infinite number of heteroclinic cycles in the network, but all 
the cycles have a common heteroclinic connection. The two networks have the same 
basic network structure as each other (see figure [LTj) but in one network, one of the 

42 



equilibria has a pair complex contracting eigenvalues while in the other network all 
eigenvalues are real; the equivariance properties of the networks are slightly different 
to accommodate this feature. 

Previous work on these and related networks pQHHHQiilllB] concen t r ated on in- 
vestigating their stability properties and understanding switching dynamics near each 
network, but did not look in detail at resonance. Here we have focussed on under- 
standing the dynamics resulting from one or more of the heteroclinic cycles in the 
network undergoing a resonance bifurcation. We have been primarily interested in 
understanding how much of the observed dynamics can be thought of as arising from 
resonance of a single cycle and how much is inherently due to the network structure. 

Our network with only real eigenvalues (Case I) contains two distinguished hete- 
roclinic cycles, one each in the subspaccs defined by j/3 = and X3 = 0. We defined 
5x (resp. 5y) to be the ratio of contracting to expanding eigenvalues seen by the cycle 
in the 2/3 = (resp. X3 = 0) subspace, and investigated the dynamics that occur for 
5x and 5y near one. When 6x or 5y passes through one, the corresponding cycle 
undergoes a resonance bifurcation and, as expected from previous work on such bi- 
furcations , a periodic orbit appears in the corresponding subspace (see 
figure [5?8| . Within each subspace, there is a transfer of stability between the hetero- 
clinic cycle and the bifurcating periodic orbit, as normally expected for resonance of 
single cycles. However, because of the network structure, none of the heteroclinic cy- 
cles can be asymptotically stable within the full phase space. This observation might 
lead one to conclude that the bifurcating periodic orbit can never be asymptotically 
stable, but we show this is not the case; the bifurcating periodic orbit may in some 
circumstances be asymptotically stable even though the cycle from which it bifurcates 
is never asymptotically stable. 

In addition to the periodic orbits that appear in the subspaces when one or other 
of the distinguished cycles goes through resonance, there may be further periodic 
orbits appearing as 5y is decreased through one, as shown in figure ISTST b) . These 
extra periodic orbits are guaranteed to exist if the quantity we called 5m, which is the 
maximum ratio of contracting to expanding eigenvalues encountered along any cycle 
in the network, is greater than one when 6y = 1. 

Resonance in the network with complex contracting eigenvalues at one equilib- 
rium (Case II) is significantly more complicated than for the case with real eigenval- 
ues. By contrast with Case I, the symmetry properties of this network do not induce 
the existence of three-dimensional subspaces in which there are distinguished hete- 
roclinic cycles. We can, however, still write down two distinguished combinations of 
eigenvalues, corresponding to two particular cycles: Sx (resp. Sy) is now the ratio of 
contracting to expanding eigenvalues seen by the orbit that approaches X (resp. Y) 
with rate determined by the contracting eigenvalue cc(0) (resp. cc(-|)) as defined in 
equations (|3.10[) and (|4.8p . We investigate the dynamics that occurs for 5x and 5y 
near one. We find that an infinite sequence of saddle-node bifurcations of periodic 
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orbits accumulates on each of the lines Sx = 1 and 5y — 1 in the (5x,5y) parameter 
space (sec figure I5.14p , and expect that there may be period doubling cascades of the 
orbits created in the saddle-node bifurcations. Note that in the Case II network, the 
quantity 5m (as defined above for the Case I network) is again always greater than the 
maximum of Sx and 5y and thus 5 m > 1 in a neighbourhood of 5x = 1 and 5y = 1. 
However, 5m niay pass through one in the region where 5x < 1 and 5y < 1. We have 
shown that the infinitely many periodic orbits created in the resonance bifurcations 
at 5x = 1 and 5y = 1 will persist so long as 5m > 1- 

In [18] , the possibility of chaotic attractors occurring in the Case II network when 
5x < 1) 5y > 1 was discussed; here we are able to estimate the location of such an 
attractor under certain conditions on the spread of orbits. In a numerical example, 
we found three co-existing chaotic attractors in the regime 5x < 1, 5y < 1. One of 
these attractors seemed to satisfy the spread condition on orbits, and its location was 
consistent with our prediction. 

Analysis of the dynamics of an isolated heteroclinic cycle with placement of the 
complex eigenvalues being analogous to the cycles in the Case II network showed (in 
section the existence of an analogous sequence of saddle-node bifurcations. We 
thus conjecture that the existence of infinitely many saddle-node bifurcations in the 
Case II example is due to the presence of the complex eigenvalues rather than arising 
from the network structure. Note that all equilibria on the isolated cycle analyzed 
in section |6] had one-dimensional unstable manifolds, and so the results from that 
example do not carry over directly to our network example, meaning we are unable 
to make a statement stronger than a conjecture at this stage. 

The bifurcations of periodic orbits we have located in our analysis appear to be es- 
sentially just those that arise from resonance bifurcations of single heteroclinic cycles, 
and provide little evidence for the effect of the network structure on the dynamics. 
However, we have restricted attention to periodic orbits that make just one circuit of 
the network before closing; it may be that orbits that make two or more circuits of the 
network (corresponding to orbits of period two or higher in the return maps) are more 
influenced by the network. One way in which the effect of the network is manifested is 
in the the role of the quantity 5m- As discussed in [18] in the context of Case II, net- 
work stability is determined by the maximum and minimum ratios of contracting and 
expanding eigenvalues experienced by any cycle in the network; the network ceases 
to be asymptotically stable when the minimum ratio (called 5 mm in [18]) decreases 
through one, and the possibility that orbits not on the stable manifold of an equilib- 
rium of the network might be attracted to the network is erased when the maximum 
ratio (called 5m here and <5 max in [18]) decreases through one. In general, neither the 
maximum nor minimum ratio is 5x or 5y but is rather some combination of eigenval- 
ues seen on different cycles. In this sense, the important combinations of eigenvalues 
for resonance of a network carry information about the network as a whole, not just 
about single cycles within the network. We note, however, that in our examples, 
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because of the geometry of the networks, the minimum ratio of eigenvalues is always 
cither 5x or Sy. 

The method of analysis we have adopted in this article is based on the standard 
procedure for construction of return maps that approximate the dynamics near the 
network, but with significant adaptations to accommodate the two-dimensional un- 
stable manifolds that occur for some equilibria; elements of the new techniques we 
have developed were first described by us in Q~8] but are extended in this article. We 
believe that similar techniques might be used for the analysis of other heteroclinic 
networks, and in particular for other networks in which all cycles have a common 
heteroclinic connection, as is the case in the two networks we considered. Analysis of 
such networks has, to date, been largely restricted to examining the dynamics near 
specific cycles in the network, but our techniques enable us to capture the dynamics 
of the whole network. The issue of extending our techniques to the study of networks 
in which cycles do not all have a common connection is left for future work. 

Finally, we note that numerical work on networks such as those considered here 
is extremely delicate. The type of analysis we have performed is, as usual, valid 
in the limit of being close to the network; we have had to look within a distance 
of 10 -60 of the network to see some of the phenomena of interest in our numerical 
examples. On the other hand, very close to our Case II network, the vast majority 
of orbits visit equilibrium X rather than Y and so it is necessary to wait for a long 
time before a typical orbit will explore the parts of the network passing near Y . 
A further complicating factor is that 8x and Sy have to be rather close to one for 
some phenomena to be observable; otherwise contraction onto or expansion away 
from the network is too rapid. Thus, while we have located a variety of phenomena 
by theoretical means, verifying the existence of all these phenomena in particular 
examples might not be straightforward. 
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